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We numerically solve 2+lD effective kinetic theory of weak coupling QCD under longitudinal 
expansion relevant for early stages of heavy-ion collisions. We find agreement with viscous hy¬ 
drodynamics and classical Yang-Mills simulations in the regimes where they are applicable. By 
choosing initial conditions that are motivated by color-glass-condensate framework we find that for 
Qa = 2GeV and as = 0.3, the system is approximately described by viscous hydrodynamics well 
before r < 1.0 fm/c. 


INTRODUCTION 

In the weak coupling picture of the pre-thermal evo¬ 
lution of heavy-ion collisions, the post collision debris 
that end up in the mid-rapidity region undergo several 
stages that are characterised by widely different physics 
and degrees of freedom. On the one hand, according to 
the saturation paradigm [1], at very early times r ^ 
where Qs is the typical energy scale right after the colli¬ 
sion, the energy is deposited in strong color helds. These 
strong helds admit a description in terms of classical 
Yang-Mills theory to leading order in ’t Hooft coupling 
A = 47rAfcOfs- Indeed, there have been several interesting 
studies of classical Yang-Mills helds under longitudinal 
expansion in the past years [2-6]. On the other hand, 
once the system has reached a local thermal equilibrium, 
or at least approximately isotropized [7], the matter in 
the mid-rapidity region (at sufficiently low pp) is de¬ 
scribed by relativistic fluid dynamics [8, 9]. There has 
been a sizeable and very successful program of numerical 
simulations of relativistic hydrodynamics. 

It is well known that classical Yang-Mills theory can 
not reach thermalization due to the Rayleigh-Jeans catas¬ 
trophe nor can it isotropize when exposed to rapid lon¬ 
gitudinal expansion as noticed by [10]. Instead classical 
evolution drives the system further away from equilib¬ 
rium, making the system less occupied (i.e. weaker helds) 
but more anisotropic, which has also been observed in 
the simulations of [4, 5]. Therefore, there is a missing 
link between the physics of saturated gluon helds and 
fluid dynamics that needs to be bridged in order to fully 
carry the predictions of saturation physics to the hydro- 
dynamic regime. This gap can be hlled with the system¬ 
atically improvable framework of effective kinetic theory 
(EKT) of [11]. 

The EKT faithfully describes to leading order in A/ 
systems where the typical occupancies of gluons are not 
nonperturbative / <C 1/A and have momentum signih- 
cantly larger than the in-medium screening scale > 
w? = X Jpf{p)/p- At weak coupling, these conditions 
are certainly fulhlled in thermal equilibrium and there¬ 


fore we can describe the system all the way to the equi¬ 
librium. While these conditions are not fulhlled at the 
very earliest times, both EKT and classical Yang-Mills 
give an equally valid leading order description for a wide 
range of large but perturbative occupancies 1 <C / 1/A 

[12, 13]. Therefore a possible strategy to simulate the 
system through all time scales is to start the simulation 
with a classical Yang-Mills simulation (for QsT ^ 1 and 
/ ~ 1/A) and subsequently pass the system to EKT at 
some arbitrary time tektQs ^ 1- Then, once EKT 
has brought the system sufficiently close to the thermal 
equilibrium, both hydrodynamics and EKT should give 
equivalent descriptions and the system can be passed 
to a hydrodynamical simulation at some arbitrary time 

'^hydro' 

The proof of principle of such a procedure was shown 
in a series of papers [13-15] in an isotropic setting in 
the absence of expansion. Here, we present first re¬ 
sults of numerical simulations of the EKT in a boost 
invariant 2-l-lD setting and demonstrate the connection 
with both classical Yang-Mills simulations and viscous 
hydrodynamics. We make a first attempt to model the 
early stages of heavy-ion collision starting from the over¬ 
occupied region all the way to the thermal equilibrium 
and find that after time QsT ~ 5.0 the time evolution 
of energy density is described by viscous hydrodynamics 
to a better than 10% accuracy for a realistic coupling 
A = 10 corresponding to ag ~ 0.3. 

The current result is, however, incomplete in two ways. 
Firstly, there have been significant steps to describe the 
melting of the strong color-fields in classical statistical 
field theory [2, 3, 6], but currently we do not have a re¬ 
liable initial state from a 3-I- ID simulation to enter into 
our EKT simulation. Therefore we initialize our system 
at QsT ^ 1 with the energy density given by 2-l-lD sim¬ 
ulation [28] and vary the parameters of the initial con¬ 
dition to quantify the ignorance of the very early time 
dynamics. Secondly, certain non-perturbative chromo- 
Weibel instabilities [16] arising from anisotropic screening 
play a parametrically leading order role during the whole 
non-equilibrium evolution [17, 18]. However, numerical 
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simulations of classical fields [6] have not been able to 
see them beyond the very early times, suggesting that 
even if parametrically leading order, their effect is numer¬ 
ically small for values of A that are phenomenologically 
interesting. Therefore, as a formally correct treatment of 
anisotropic screening is rather complicated [19-22] and 
not a fully solved problem, we in this first work will treat 
the screening in an isotropic way. Therefore, our results 
are correct to leading order in A for systems that are close 
to isotropy (late times) and for large anisotropies our re¬ 
sults have leading logarithmic accuracy apart from the 
instabilities. 


METHODOLOGY 


The EKT of [11] is defined though the effective Boltz¬ 
mann equation for the color and spin averaged distri¬ 
bution function of gluons, and to leading order in A it 
contains effective 2 0 2 scattering and 1 o 2 splitting 
terms. 

= Cl-H-2[/p] + C2o2[/p] + Cexp[/p]- (1) 

We restrict ourselves to azimuthally symmetric distri¬ 
butions but allow anisotropy in the z-direction so that 
it is enough to specify /p = fxp,p with Xp = z • p. 
The effect of longitudinal expansion is encapsulated in 
Cexp[/](p) = [23]. 

The 2 o 2 effective scattering term reads 

C2„d/Kp) = / dTMMf 

X (/p /k <7p' 5k' ~ fp' /k' 5p 5k) (2) 

X [(5(p - p) -b S{p - k) - 5(p - p') - S{p - k')], 


where v = 2d a and gp = 1 + fp. dTps is the integral 
measure over the phase space of 2 o 2 processes, singling 
out the z-direction 


dTps = 


211^7 


dq / duj 
J —q . 
/*! p 2 tt 


/ dp dk 

{q-uj)/2 J{q+ui)/2 


dx„ 


dcjdpqdcjij^q. 


(3) 


In these coordinates the angles of incoming and outgoing 
momenta read 

X|p} sin cos \/^ d^p'^qXq . (4) 


elation holds 

for Xk, 

Xpi, and 

Xk' with cos cj)p'q = 

jr, COS (pk'q 

= cos (j)kq, 
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and t = uj'^ — . The effective matrix element for most 

kinematics is the ordinary vacuum tree-level element 


\M\^ 

WX'^dA 


9 {s-tr ju-s)^ 

4 m2 + ^2 


S2 ) 


(7) 


For small t ~ or m ~ the tree-level element 
has an IR divergence that is regulated by the physics 
of screening. In the formulation of [11], the screening is 
treated by replacing the naive tree-level matrix element 
by the retarded Hard Thermal Loop (HTL) resummed 
expression in the soft kinematic region. The resummed 
propagator in the anisotropic case, however, contains a 
non-integrable singularity, signalling the presence of an 
instability in the system. Therefore, in the anisotropic 
system this procedure is not fully satisfactory. Here, in 
this work we will restrict ourselves to isotropic screening 
and use the prescription derived in [13], replacing in the 
denominator (similarly for u) 

t(g2-b2^gm^), (8) 


with ^0 = -x/S, which is leading order accurate in the 

case of isotropic distribution. 

The effective 10 2 collinear splitting term reads 

(0Tr]3 1 /-oo i.p/2 

Ci«2[/](p) = / dp / dk'[47T^{p;p',k')] 

47rp2 p Jq 

X {fxx,p9x^,p'gx^,k' - 9xp,pfxp,p'fxp,k') 

X [d{p - p) - 6{p - p') - S{p - k')] , (9) 


where p' = p — k\ and the effective splitting rate reads 


7(p;p',fc') 


p'^ + p''^ + k"^ dAX 

p3p/3k/3 2(27r)^ 


d^h 


2h • ReF. 


Again, consistently assuming only isotropic screening, 
the function F is given by the solution to the linear inte¬ 
gral equation 


2h = iSE{h)F{h) + 


d^q 

(^)2L' 


-4(q) 


( 10 ) 


X (3F(h) — F(h — pq) — F(h — fcq) — F(h -f p'q)) 


= 2^ ! + /p)’ - 

mf Ip -b /2pk'p'. We solve this equation using an effi¬ 
cient numerical method introduced in [24]. 

For the numerical solution of the FKT, we discretize 
'nxj,,p = 47rp^/(27r)^/xp,p on a 2D grid and Monte Carlo 
estimate the integrals of Fqs. (3,9). We use a logarith¬ 
mically spaced grid for both p and Xp with at least 250 
grid points in angular direction and at least 100 in the 
p-direction. We have varied the number of grid points 
to verify that the results are insensitive to the number of 
grid points. Our algorithm, the discrete-p method of [13], 
conserves energy exactly and has exactly correct particle 
number violation for 1 o 2 term. 
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FIG. 1: Trajectories of runs with different initial conditions 
^ = 4 (Solid lines) and ^ = 10 (dahsed lines) and varying 
coupling A in a plane of mean occupancy (weighted by the 
energy of particles) and anisotropy. The A = 0 line corre¬ 
sponds to classical field approximation. The violet dots refer 
to the times in Fig. 2. The simulations at finite coupling reach 
thermal equilibrium located at points indicated by the black 
crosses. 


RESULTS 


We shall now apply EKT to simulate the prethermal 
evolution of the expanding fireball created in a heavy- 
ion collision. In saturation framework, the initial condi¬ 
tion is typically described in terms of “gluon liberation 
coefficient” c and mean transverse momentum (jpt) jQs 
[26, 27]. The gluon liberation coefficient is proportionaf 
to the total gluon multiplicity per unit rapidity 


2dAT 


(Pp _ dNjnit.g dAQl 

(27r)3 d’^x±dy ttA 


( 11 ) 


after the classical fields have decohered and can be de¬ 
scribed in terms of quasi-particles. Lappi [28] finds in 
JIMWLK evolved MV model values relevant for heavy- 
ion collisions relevant for LHC of roughly {px) ~ 1.8Qs 
and c ~ 1.25 extracted at time QgT = 12 from a 2D 
classical Yang-Mills simulation. By construction the dis¬ 
tribution then has {pz) = 0. But it is has been noted [25] 
that certain plasma instabilities will broaden the distri¬ 
bution in pz in a time scale Qt ~ l/log^(A“^). There¬ 
fore, as a rough estimate of the initial condition we in¬ 
stead take somewhat arbitrarily our initial condition at 
the time Qt = 1 to be 

f{Pz,Pt) = jAfo{pz^/{pT),P±/ (pt)), (12) 

fo{p.,PA) = ^=L=e-2(pi+P^)/3, (13) 

Vp±+pi 


Figure 1 displays a set of trajectories from simula¬ 
tions with varying A and ^ = 4,10 on a plane of mean 
occupancy (weighted by the energy of particles) and 
anisotropy measured by the ratio of the transverse and 
longitudinal pressures Pt/Pl- The line with A = 0 cor¬ 
responds to the classical field limit A —?► 0 with fixed A/, 
which is obtained in EKT by including only the highest 
power of /’s in Eqs.(2,9). The classical field theory can 
not thermalize and indeed instead it flows to a stationary 
scaling solution. By performing classical Yang-Mills sim¬ 
ulations Berges et. al have established that the scaling 
solution can described by a scaling form of the distribu¬ 
tion function [4], 

f{Pz,P±,T) = {QsT)~'^^^fs{{QsT)^^^Pz,P±), (14) 

where fs is approximately constant as a function of time. 
This behavior is demonstrated in Figure 2 where we plot 
a section of rescaled distribution function fs measured 
at various times as a function for pz = {QsTy^^Pz at 
fixed p± following Berges et ah. Our results corroborate 
that such a scaling solution exists at late times within the 
classical approximation and we observe that the scaling 
regime is reached after a time QsT ^ 15. 

Moving on to the finite but small coupling A = 1,0.5, 
we see qualitative agreement with the parametric picture 
of bottom-up thermalization of [10]: there are three dis¬ 
tinct stages of evolution visible. In the first stage the 
classical evolution drives the system more anisotrpic and 
less occupied. Once the occupancies reach / ~ 1, there 
is a qualitative change in the dynamics of the system as 
Bose enhancement is lost. This has the effect that an- 
siotropy freezes but the system still continues to get more 
dilute. Only in the last stage which is characterized by 
a radiational break up of the particles at the scale Qs, 
does the trajectory turn back and reach thermal equi¬ 
librium, denoted by the black crosses in the Figure 1. 
For larger values of coupling A = 5.0,10, these features 
become however less pronounced and the system takes 
more straight trajectory towards equilibrium. It should 
be noted that for these values of A the assumption of 
p ^ m is not satisfied and large NLO corrections are to 
be expected. 


Approach to hydrodynamics 

We expect that late time evolution should be described 
by relativistic hydrodynamics. Under flow with transla¬ 
tional invariance along transverse directions and boost 
invariance, the hydrodynamical relations read to second 
order in gradients [29] 


choosing A such that comoving energy density re = 
{pT)dN/d‘^xdy is fixed. We then vary ^ = 4,10 to quan¬ 
tify our ignorance of the initial nonperturbative dynam- 
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FIG. 2: Sections of scaled distribution fs{Pz,P±) = 
iQsr)^^^fiPz{Q3'r)~^^^,p±) at p± = l.SQs in classical ap¬ 
proximation at vastly different times. The good overlap of the 
curves indicates that system has reached the classical scaling 
solution of Eq. (14). In contrast, QsT = 5 has not yet reached 
the scaling solution. 


with longitudinal and transverse pressure Pp = — $ 

and Pt = First order hydrodynamics corre¬ 

sponds to setting $ = in Eq. (15). At weak coupling, 
the transport coefficients rj, rn and Ai are known [30, 31] 
leaving zero free parameters to fit. besides a time when 
the hydrodynamics is initialized. We initialize the 1st 
order hydrodynamics at the latest time we have in our 
simulation and integrate Eq. (15) backwards in time. For 
2nd order hydrodynamics integrating backwards is highly 
unstable and we initialize the energy density at some ar¬ 
bitrary earlier time and integrate forwards in time. 

In Figure 3 we examine the validity of the hydrody- 
namical expansion at small A = 1 and at realistic in¬ 
termediate A = 10 (a^, ~ 0.3) values of coupling. In 
both cases we see that the evolution of the components 
of the energy momentum tensor asymptotes to their hy- 
drodynamical values. In case of A = 1, the hydrody- 
namical behaviour is reached only at a rather late time 
QsT ~ 2000. We have checked that including 2nd order 
terms before this time does not make the convergence 
significantly better; before this time the evolution differs 
qualitatively from the hydrodynamical prediction. How¬ 
ever, rather remarkably, for A = 10 even 1st order hy¬ 
drodynamics gives a very good description of the data 
all the way to very early times QsT ^ 10 (corresponding 
to T ~lfm/c for Qs = 2.0GeV) where the ratio of the 
pressures is still as large as Pt/Pl ~ 5. In addition, in¬ 
cluding the second order terms significantly improves the 
convergence. Indeed, we find that initializing the 2nd or¬ 
der hydrodynamics at QbT = 1 leads to only 10% relative 
error in the energy density at late times. 


DISCUSSIONS AND CONCLUSION 

The parametric estimate of Baier et al. [10] for the 
time when the hydrodynamic behaviour sets is QsT ^ 




FIG. 3: Longitudinal pressure Pl, energy density e, and trans¬ 
verse pressure Pt from a simulation with ^ = 10.0 and A = 1 
(top panel) and A = 10 (bottom panel). The components of 
Tj^ have been scaled by so that the ideal hydrodynam¬ 
ics corresponds to horizontal lines. The scale on x-axis with 
[fm/c] corresponds to Qs = 2.0GeV« 10/fm. The simulations 
with 5 = 4 are also displayed with thin dotted lines. 



A — 1/3, . \" 4 / 3 ,„ , 

A, (Ti/s) (Qx) 

FIG. 4: Scaling of the approach to 1st order hydrodynamics 
at various values of A. The dashed lines correspond to 5 = 4 
whereas the solid ones have 5 = 10. 


A“13/^. This arises from equating the age of the system 
T with the time scale tq it takes to affect appreciably 
the scale Qs in a thermal bath whose temperature de¬ 
pends on this time T{t) ^ \~^^'^Qs[Qst)~^^'^ according 
to conservation of comoving energy density. In [10] it 
was assumed that the rate for affecting the scale Qs is 
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LPM suppressed [32], giving tq ^ {QsITY/"^/\^T. A 
selfconsistent solution of these equations gives the afore¬ 
mentioned estimate of [10]. Arnold and Lenaghan [33] 
noted that as the scattering rate in thermal plasma is 
tq ~ 1/A^T, there can be no process that would make 
the estimate faster than QsT ^ We have exam¬ 

ined the validity of these scaling estimates by plotting the 
difference of the energy density obtained from the simu¬ 
lation and the 1 st order hydrodynamic estimate and find 
that both of these estimates describe the data poorly. 
However, if we assume that the approach is governed by 
the hydrodynamical relaxation time tq ~ rn ~ rj/sT 
we get an estimate r ~ /Qs. Figure 4 dis¬ 

plays the deviation of the hydrodynamics as a function 
of rescaled time. In particular for intermediate couplings 
A = 5,10 the overlap of the different curves indicates 
correct scaling. Note that this estimate is parametri¬ 
cally the same as the estimate of Arnold and Lenaghan 
as parametrically ij/s ^ However, there are large 

corrections beyond the parametric estimate in 77 /s due 
to which it is important to use the full numerical result 
instead of the simpler parametric estimate. We however 
believe (in the absence of plasma instabilities) that the 
correct scaling at sufficiently small A is that of [10]. This 
estimate is however based on a large scale separation of 
T(Thydro)/Qs A^/®. Numerically this ratio is not very 
large even in our simulation with small A = 1.0, where 
the ratio is T/Qs Ri 0.15 (at QsT = 1800) as can be in¬ 
ferred from Figure 3. Therefore we believe that the scal¬ 
ing predicted by [ 10 ] sets in only at significantly smaller 
couplings. 

In conclusion, we have shown how a far from equi¬ 
librium overoccupied configuration of gluons reaches hy¬ 
drodynamical flow under longitudinal expansion in a 
weak coupling setting that is systematically improvable. 
It has been demonstrated by Chesler and Yaffe within 
AdS/CFT framework that at large values of’t Hooft cou¬ 
pling, hydrodynamics is surprisingly robust even in the 
presence of large anisotropy [34]. The main deliverable 
of this paper is to show that this robustness is present 
also in the weak coupling picture extrapolated to inter¬ 
mediate couplings relevant for heavy-ion collisions. 
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